I want to try and probe a question that has been raised since Las Vegas, can a state regulate firearms on it’s own, given the free trade between jurisdictions. Since there is not an open electronic federal database for firearm ownership and transactions, one ipportant source of information is the Bureau of Alcohol, Tobacco and Firearms, (ATF). They publish the trace of firearms that were recovered every year, and when possible trace the state where the firearm originated. This creates a matrix that is similar to what in economics is called Direction of Trade.

In this matrix there is the rows depict the source and the columns depict the destination. This lets one get an idea where firearms that are confiscated by the ATF orginated from. From this we can also infer which states are net importers of firearms and which states are net exporters.

I will explore this matrix in an attempt to better understand if firearms are more likely to flow between geographically adjacent states. In the end I will get to a shiny app that ties everything together, for those who want to stop here…

Here is the a short gif showing the app

Below is the script to run the app from R


pkgs <- c('reshape2','geojson','readxl','ggplot2',
'leaflet','httr','rgeolocate','shiny','sp','dplyr')

check <- sapply(pkgs,require,warn.conflicts = TRUE,character.only = TRUE)

if(any(!check)){
  pkgs.missing <- pkgs[!check]
  install.packages(pkgs.missing)
  check <- sapply(pkgs.missing,require,warn.conflicts = TRUE,character.only = TRUE)
  }

shiny::runGitHub('yonicd/gunflow')

So we start by loading all the pacakges we need

pkgs <- c('reshape2','geojson','readxl','ggplot2','heatmaply',
'leaflet','httr','rgeolocate','shiny','sp','dplyr')

sapply(pkgs,require,character.only = TRUE)

Read in the data

f0 <- tempfile()
download.file("https://raw.githubusercontent.com/rstudio/leaflet/gh-pages/json/us-states.geojson",destfile = f0)
states <- geojsonio::geojson_read(f0, what = "sp")

f1 <- tempfile()
download.file('https://www.atf.gov/docs/undefined/sourcerecoverybystatecy2016xlsx/download',destfile = f1)

gun_mat <- readxl::read_xlsx(f1,col_names = TRUE,range = 'B2:BD56')%>%
data.frame()

A natural place to inspect a matrix, a heatmap. We clean up the data.frame to set it up to use the heatmaply package.

g <- gun_mat[,-1]
row.names(g) <- gun_mat[,1]

mg <- as.matrix(g)

#remove diagonal (source of firearm in the state itself)
diag(mg) <- NA

mg <- mg[!grepl('^GUAM|^US VIRGIN',rownames(mg)),!grepl('^GUAM|^US VIRGIN',colnames(mg))]

heatmaply::heatmaply(heatmaply::percentize(mg),
          main='Firearms Sourced and Recovered in the United States and Territories 2016',
          xlab = 'Source',
          ylab = 'Recovered',
          colors = rev(heatmaply::RdYlBu(256)),
          k_row = 5,
          k_col = 5,
          fontsize_row = 6,
          fontsize_col = 6)

We can see that geographically adjacent states are placed in clusters (west coast, east coast, et. al.), and they exhibit high levels of similarity. But this still keeps the picture pretty complex.

We move on to geographic based visualizations, using the leaflet package. Using the maps we can set a state as the state of interest and check what states are recieving firearms from it or what states are supplying it with firearms. Below is an example of such a map where Ohio is the state of interest. This map shows the outflow of firearms which is transformed as a percentage of total firearms (excl Ohio itself).

Here we see that firearms that originated in Ohio found there way to all their adjecent states (MI, IN, KY, PA, WV), but also non-adjecent ones (NY, NC, FL, CA, IL,TX).


# Function that capitalizes first letter in each word
capitalize=function(x){
  gsub("(^|[[:space:]])([[:alpha:]])", "\\1\\U\\2", x, perl=TRUE)
}

#reshape the matrix into long format and organize names
gun_mat <- gun_mat%>%
  reshape2::melt(.,'X__1')


names(gun_mat) <- c('from','to','value')

gun_mat <- gun_mat %>%
  dplyr::filter(!(grepl('^GUAM|^US VIRGIN',to)|grepl('^GUAM|^US VIRGIN',from)))%>%
  dplyr::mutate_at(.vars = vars(to,from),.funs=funs(capitalize(tolower(.))))

#Set state and flow direction
thisstate <- 'Ohio'
type <- 'Outflow'

#Calculate percent changes needed
gun_mat1 <- switch(type,
                   Inflow={
                     gun_mat%>%
                       dplyr::group_by(to)%>%
                       dplyr::mutate(value=ifelse(to==from,NA,value),pct=100*value/sum(value,na.rm = TRUE))%>%
                       dplyr::filter(to==thisstate)%>%
                       dplyr::rename(state=from)       
                   },
                   Outflow={
                     gun_mat%>%
                       dplyr::group_by(from)%>%
                       dplyr::mutate(value=ifelse(to==from,NA,value),pct=100*value/sum(value,na.rm = TRUE))%>%
                       dplyr::filter(from==thisstate)%>%
                       dplyr::rename(state=to)
                   })

#merge back into the states object
mydata <- states@data   
mydata <- mydata%>%
  rename(state=name)%>%
  mutate(state=as.character(state))%>%
  left_join(gun_mat1%>%ungroup%>%select(state,value,pct),by='state')

states@data$pct <- mydata$pct
states@data$level <- mydata$value
states@data$density <- NULL

#set pallete
pal <- colorNumeric(
  palette = "RdYlBu",
  domain = states@data$pct,na.color = 'black',reverse = TRUE)

#Create hover labels
labels <- switch (type,
                  Inflow={
                    sprintf(
                      "Of the %s Out of State Firearms Recovered in <strong>%s</strong><br/> %g%% of them originating from <strong>%s</strong>",
                      sum(states@data$level,na.rm = TRUE),
                      thisstate,
                      round(states@data$pct,2),
                      states$name
                    )
                  },
                  Outflow={
                    sprintf(
                      'Of the %s Out of State Firearms Originating from <strong>%s</strong><br/> %g%% were Recovered in <strong>%s</strong>',
                      sum(states@data$level,na.rm = TRUE),
                      thisstate,
                      round(states@data$pct,2),
                      states$name
                    ) 
                  }
)%>% lapply(htmltools::HTML)

#build widget

Sys.setenv(MAPBOX_ACCESS_TOKEN='AN ACCESS TOKEN FROM MAPBOX')

m <- leaflet(states) %>%
  setView(-96, 37.8, 4) %>%
  addProviderTiles("MapBox", options = providerTileOptions(
    id = "mapbox.light",
    accessToken = Sys.getenv('MAPBOX_ACCESS_TOKEN')))%>% 
  addPolygons(
  fillColor = ~pal(pct),
  weight = 2,
  smoothFactor = 0.2,
  stroke=FALSE,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 1,
    bringToFront = TRUE),
  label = labels,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto"))%>% 
  addLegend(pal = pal, values = ~pct, opacity = 0.7, title = 'Percent',
            position = "bottomleft",na.label = 'Selected State')

The final plot that is created is the net flow for each state (the sum of inflow minus the sum of outflow). This can be a problem for big states like California, they will look like they have a high net inflow, just because of their size. To take this into account the flow per 100 firearms is calculated for each state and then the inflow rate is subtracted from the outflow rate.


capitalize=function(x){
  gsub("(^|[[:space:]])([[:alpha:]])", "\\1\\U\\2", x, perl=TRUE)
}


calc <- function(side){
  total <- gun_mat%>%
    group_by_(side)%>%
    summarise(total_sum=sum(value))%>%
    rename_('state'=side)
  
  nonstate <- gun_mat%>%
    filter(to!=from)%>%
    group_by_(side)%>%
    summarise(state_sum=sum(value))%>%
    rename_('state'=side)
  
  ret <- nonstate%>%
    left_join(total,by='state')%>%
    mutate(ratio=100*state_sum/total_sum)
  
  names(ret)[-1] <- paste(names(ret)[-1],side,sep = '_')
  
  ret
}

net_flow <- calc(side = 'from')%>%
  left_join(calc(side = 'to'),by='state')%>%
  mutate(net=state_sum_from-state_sum_to,
         ratio_net=ratio_from-ratio_to)%>%
  arrange(desc(ratio_net))

net_flow$state <- factor(net_flow$state,levels = net_flow$state)

ggplot2::ggplot(net_flow,ggplot2::aes(x=state,y=ratio_net,
                                      fill=cut(ratio_net,breaks = 10,include.lowest = TRUE)))+
  ggplot2::geom_bar(stat='identity')+
  scale_fill_brewer(palette = "RdYlBu",direction = -1,name=NULL)+
  theme_bw()+
  labs(title='Net Firearm Flow per 100 Firearms Between States',
       subtitle='High is Net Exporter, Low is Net Importer',
       y='Net Ratio per 100 Firearms',x='State')+
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle=90),legend.position = 'bottom')

Last, we take all those static parts and combine them into a shiny app! Adding controls for firearm flow and attaching the hover to the set flow so arrows to show what state is selected and what state is being hovered over.

Full script of the shinyapp can be found in the github repository.